Thermodynamic first order transition and inverse freezing in a 3D spin-glass. 
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We present a numerical study of the random Blume-Capel model in three dimension. The phase 
diagram is characterized by spin-glass/paramagnet phase transitions both of first and second order 
in the thermodynamic sense. Numerical simulations are performed using the Exchange-Monte Carlo 
algorithm, providing clear evidence for inverse freezing. The main features at criticality and in the 
phase coexistence region are investigated. We are not privy to other 3D short-range systems with 
quenched disorder undergoing inverse freezing. 
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Introduction. We aim to investigate the phenomenon 
of reversible inverse transition (IT), occurring between a 
solid and a liquid in the inverse order relation relatively 
to standard transitions. The case of "ordering in disor- 
der" , occurring in a crystal solid that liquefies on cooling, 
is generally termed inverse melting. If the solid is amor- 
phous the IT is termed inverse freezing (IF). 

IT's are observed in different materials. The first 
examples were the low temperature liquid and crystal 
phases of helium isotopes He 3 and He 4 [l|. A more re- 
cent and complex material is methyl-cellulose solution in 
water, undergoing a reversible inverse sol-gel transition 
0. Other examples are found in poly(4-methylpentene- 
1) - P4MP1 - at high pressure [3], in solutions of a- 
cyclodextrine (cvCD) and 4-methypyridine (4MP) in wa- 
ter Q , in ferromagnetic systems of gold nanoparticles [|[ 
and for the magnetic flux lines in a high temperature su- 
perconductor [6|. In mentioning these cases, we stick to a 
definition of IT as the one hypothesized by Tammann Q 
a century ago: a reversible transition in temperature at 
fixed pressure [|[ whose low T phase is an isotropic fluid 
(or paramagnet - PM - for magnetic systems). IT is not 
an exact synonim of reentrance. Indeed, a reentrance can 
be absent, as for cvCD [1] or methyl-cellulose [1] solutions, 
where no high temperature fluid phase is detected. More- 
over, not all reentrances are signatures of an IT, as for 
those phases with different kind of symmetry separated 
by reentrant isobaric transition lines in temperature, cf., 
e.g., Refs. [9j, in which, however, no melting to a com- 
pletely disordered isotropic phase is present. 

A thourough explanation of the fundamental mecha- 
nisms leading to the IT would require a microscopic anal- 
ysis of the single components behavior and their mutual 
interactions as temperature changes accross the critical 
point. Due to the complexity of the structure of poly- 
mers and macromolecules acting in such transformations 
a clear-cut picture of the state of single components is 
seldom available. For the case of methyl-cellulose 
where methyl groups (MGs) are distributed randomly 
and heterogeneously along the polymer chain [Toj , Haque 
and Morris [111 ] proposed that chains exist in solution as 
folded hydrophilic bundles in which hydrophobic MGs 



are packed. As T is raised, bundles unfold, exposing 
MGs to water molecules and causing a large increase in 
volume and the formation of hydrophobic links eventu- 
ally leading to a gel condensation. The polymers in the 
folded state are poorly interacting but also yield a smaller 
cntropic contribution than the unfolded ones. A similar 
behavior has been recently modeled in colloidal systems 
by the so-called Janus particles 

Under the assumption that this is at least one of 
the fundamental mechanisms underlying IT, we model 
it approximating the folded/unfolded conformation by 
bosonic spins: s = representing inactive state, s ^ 
interacting ones, cf. Refs. [l3lll4]|. To represent the ran- 
domness on the position of the " interaction carrying" el- 
ements (e.g., MGs) we will introduce quenched disorder. 
Theoretical modeling for IT's mainly consists in heuristic 
reproductions of the phenomenon [l3l. [l5j. In particular, 
IF, has been recently observed in spin-glass (SG) mean- 
field (MF) Blume-Emery-Griffiths-Capel (BEGC) mod- 
els with spin— 1 variables 0, We will focus on the 
random Blume-Capel (BC) model [llj], whose MF solu- 
tion [l9| predicts a phase diagram with both a SG/PM 
second order and a first order phase transition (FOPT), 
i.e., displaying latent heat and phase coexistence |20| . 
The latter is characterized by the phenomenon of IT 
[l|| [2l[ . This is at contrast with the behavior of the orig- 
inal ordered BC models, in which no IT was observed 
[H-HI and with the behavior of a 3D BC model with 
quenched disorder on a hierarchical lattice (25[, yielding 
no IF, nor first order SG/PM transition. Eventually, to 
the best of our knowledge, the only claim of the existence 
of a FOPT in D = 3 systems in presence of quen ched dis- 
order, has been made for the 4- Potts glass [26|. In that 
case, though, randomness tends to strongly smoothen the 
transition into a second order one. 

Motivated by the above considerations we have, thus, 
numerically studied the existence of IF in a 3D random 
BC model with nearest-neighbor (nn) interactions. 

The model. We consider the following Hamiltonian 
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where (ij) indicate ordered couples of nn sites, and s, = 
— 1,0, +1 are spin-1 variables lying on a cubic lattice of 
size N — L 3 with periodic boundary condition. Crystal 
field D is a chemical potential for the magnetically active 
sites. Random couplings Jij are independent identically 
distributed as P(J l] ) = 1/2 5 (Jij - 1) + 1/2 5 (J^ + lj. 

We simulate two real replicas {s\^} and {s^ } of the 
system and define the overlap, order parameter of SG 
transition, as = l/iV^fsf sf ), where (. . .) is the 
thermal average. If a FOPT occurs, the order param- 
eter characterizing the transition is the density of mag- 
netically active (|sj| = 1) sites: — l/^X)i( s ?) or ' 
since we deal with finite size (FS) systems, its distribu- 
tion Pn,j(p)- The values of the parameters depend on 
the particular realization of disorder ({Jij})- Such de- 
pendence is self-averaging for t he densit y probability dis- 
tribution (Pnj(p) ~ Pn(p) = Pn,j(p) for N > 1), but 
not for the overlap distribution. We denote by ... the 
average over quenched disorder. 

Finite Size Scaling (FSS) for Continuous Transi- 
tions. In order to infer the details of the critical behavior 
from numerical simulations of FS systems, a fundamen- 
tal quantity is the four-point correlation function, i.e., 
the correlation between local overlaps qt = s^ s[ 2 ^: 



Ci(r) = Jf^2(<li Qi+r) 



(2) 



The information contained in C4 can be exploited to iden- 
tify the existence of a second order phase transition for 
FS systems, e.g., looking at a FS correlation length-like 
scaling function defined, on a 3D lattice, as (27j : 



e c = 



1 



( &(0) 

4sin 2 /ci/2 1 C 4 (fci) 



(3) 



where Cn(k) is the Fourier transform of C4(r), k\ = \k 1 \ 1 
and fc x = (^,0,0) is the minimum wave-vector. In the 
thermodynamic limit, a second order transition is charac- 
terized by a diverging correlation lenght, at critical tem- 
perature T c , whose FSS behavior is the same as in Eq. §5§ 
[28l [29j . Another relevant observable is the SG suscep- 
tibility X SG = N{f) = 64(0), diverging at the PM/SG 
transition as N — >• 00. Because of FS, though, £ c and 
Xsg cannot diverge in numerical simulations. Around 
the critical region, however, scale invariance survives. In 
fact, we can define a FS "critical" temperature T C L as the 
temperature at which the above mentioned observables 
do not depend on the size. In this scaling region we have 

UL = Utc/L) = i(L^(T - T c )) (4) 
XSG^- 2 = jfa/L) = x(lM»(T - T c )) (5) 

The critical temperature can, then, be estimated by FSS 
of T C L in the L ->• 00 limit. 

In order to estimate the critical exponents we use the 
FSS quotient method [28| , based on the observation that 



at Tjf the correlation lenghts of different linear sizes L 
and sL (in L units) are equal: s£ c (T c L , L) = £ C (T C L , sL). 
For an observable A diverging as t XA (t = T/T c — 1) this 
implies: 



A(T^sL) 
A(T£,L) 



0(L- W ). 



(6) 



For a SG we can obtain the exponents 77 by means of 
the FSS of the quotients of dp£ c and XSG, scaling, respec- 



tively, with exponents xq„ 



1 



v and x XSG = (2 



■T))v. 



Characterization of First Order Transitions. The 
Clausius-Clapeyron equation for our model, where D 
plays the role of a pressure, reads }16| 



dD _ s PM - ssg _ As 
dT ppm - Psg Ap 



(7) 



When the system undergoes a FOPT, a discontinuous 
jump in p (and, thus, in q) occurs. At finite N, Pn(p) 
displays two peaks in the coexistence region correspond- 
ing to PM (/?pm) and SG (psg) phases. The FS transi- 
tion line P C (N, T) can be evaluated as the locus of points 
where the two phases are equiprobable, i.e., the areas of 
the two peaks are equal [30| : 



dpP N (p) 



dpP N (p), 



(8) 



with po £ [ppm , Psg] such that Pn(po) = (or minimal 
for small N next to the tricritical point). 

Exchange Monte Carlo in T and P. We simulated 
the equilibrium dynamics of our model using the parallel 
tempering (PT) algorithm, replicating the system both 
on different parallel T and P. For the PT in T, the swap 
probability of two copies between T — T and T + AT 
is: P su , ap (A/3) = min [1, exp{A/3AH}]. Between P and 
P + AP it reads P swap (AP) = min [1, exp{/3ADAp}]. 
We used the latter to identify the reentrance of the tran- 
sition line in the T, P phase diagram, cf. Fig. [TJ 

We studied 3D systems with PT in T at P = 
0,1,1.75,2,2.05,2.11, and in P at T = 0.2,0.3,0.4,0.5. 
At all P we simulated from 33 to 40 replicated thermal 
baths Nt at linear size L = 6,8,10,12 (number of disor- 
dered sample: Nj = 2000). For P = 0, 1, 1.75, 2 we sim- 
ulated N T e [20 : 33] at L = 16,20 (Nj e [900 : 1500]) 
and N T G [17 : 22] at L = 24 (Nj G [500 : 1000]). For 
the PT cycles in P, No G [21 : 37], parallel replicas 
at different P were simulated, of L = 6, 8, 10, 12 and 15 
(Nj = 1000). In the latter case varying AD were used, 
larger in the pure phases and progressively smaller ap- 
proaching the coexistence region. The number of Monte 
Carlo (MC) steps varies from 2 15 to 2 21 according to L 
and to the lowest values of T, P reached. Thcrmalization 
has been checked by looking at: (i) the symmetry of the 
overlap distributions Pn,j(q), (h) the t-log behavior of 
the energy (when at least the last two points coincide), 
(iii) the lack of variation of each considered observable 
( e -g-j £,iXsg) on logarithmic time-windows. 
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FIG. 1: Phase diagram in D, T: second order transition and 
an inverted FOTP occur. In the latter case also the spinodal 
lines are reported (dashed). Bottom inset: latent heat \ As\/T 
along the first order line. Top inset: detail of IF region, inter- 
polation of transition line D c (oo,T) (dotted), spinodal lines 
(dashed). The error bars are the FSS of the minimal inter- 
val in T and D at each L needed to identify the crossings in 
£ C /L curves (for continuous transitions) or compare the areas 
under P N (p) for FOPT. 



Numerical results. In Fig. f2]we present the T-behavior 
of Z/L for D = 0,1,2,2.11. From the FSS analysis of 
their crossing points we can determine the critical tem- 
perature and, applying the quotient method (cf. Eq. © 
with s = 2), we obtain estimates for the critical expo- 
nents: v = 2.44(6), 77 = -0.34(2) at D = 0, v = 2.4(2), 
r] = -0.31(2) at D = 1, v = 2.1(2), j] = -0.27(2) at D = 
1.75. The system appears to be in the same universality 
class of the Edwards- Anderson model (corresponding to 
the D = -00 limit of our model) [H,[3l|. For D = 2, near 
the tricritical point (0.5 < T < 0.54,2.05 < D < 2.11), 
the quotient method does not yield reliable estimates be- 
cause of a crossover in the scaling functions in the range 
of probed sizes (L = 6 to 24). This comes about because 
at the tricritical point the coefficient of the fourth order 
term in the SG free energy action goes to zero and the 
sixth order term becomes relevant for the critical behav- 
lil, a typical behavior of BEGC-like systems [321 ]. 



lor 

At D — 2.11 no evidence is found for a second or- 
der phase transition whereas a FOPT is observed at 
T = 0.51(1). This is also the reason why, as shown in 
the bottom-right panel of Fig. [2] £ C /L decays in the co- 
existence region, for T < 0.51. 

The FOPT is determined by looking at Pl(p) as two 
peaks appear, and computing the (-D, T) points satisfy- 
ing Eq. JS} . A pure phase corresponds to a single-peaked 
distribution, whose peak is narrower for larger systems. 
In Fig. [3] we show the behavior of the Pis(yo) through the 
FOPT in D at T = 0.4 [33]. A detail of the FSS limit of 
the IT line D C (T) and of the spinodal lines is plotted in 
the bottom inset of Fig. [T] The spinodal lines at given 
L are estimated by looking at the D values at which a 
secondary peak first arises. 
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FIG. 2: Scaling funtions £ C /L vs. T for different val- 
ues of the chemical potential D. For D = 0, 1, 2 (L = 
6, 8, 10, 12, 16, 20, 24) a continuous phase transition is found in 
the region of scale invariance. At D = 2.11 (L = 6, 8, 10, 12) 
no crossing is observed and at low T £ C /L 0. 
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FIG. 3: Density distribution Pl(p), L — 15, across the co- 
existence region at T = 0.4: two peaks develop at ppm and 
Psg- As D increases the thermodynamically relevant phase 
(lowest free energy) passes from SG to PM in a first order 
phase transition. The dominant phase correpsonds to the one 
with larger probability, i.e., larger intergral of the peak. As 
the peak at psa vanishes the system is in a purely PM phase. 
Inset: Pis(p) on t/-Log scale. 



Using Eq. 0, from the knowledge of Ap and the es- 
timate of dD I dT by numerical interpolation we compute 
the latent heat \As\/T — (ssg — s pm)/T employed in the 
transition, cf. bottom inset of Fig. [TJ as T increases the 
PM acquires latent heat to vitrify; SG entropy is higher 
than spm and the "frozen" phase is found at a higher T 
than the fluid one. The IF takes place between a SG of 
high density to an almost empty PM (e.g., at T=0.4, 
in the coexistence region D G [2.046(2) : 2.092(5)], 
Psg — 0.52 and ppm — 0.03). The few active sites 
do not interact with each other but only with inactive 
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neighbors and this induces zero magnetization and over- 
lap. The corresponding PM phase at high T has, instead, 
higher density (e.g., p PM {D = 2,T = 0.6) = 0.4157(2), 
p PM (D = 2.11, T = 0.6) = 0.596(2)) and the param- 
agnetic behavior is brought about by the lack of both 
magnetic order (zero magnetization) and blocked spin 
configurations (zero overlap). 

Conclusions. We focused on a spin-1 SG model on a 
3D cubic lattice, whose constituent features try to cap- 
ture at least one supposed mechanism underlying IT: the 
raise of inactive components at low T. We provide nu- 
merical evidence for an equilibrium inverse freezing phe- 
nomenon: at given values of an external field, heating up 
a paramagnet this is transformed into a SG. The whole 
phase diagram in temperature T and crystal field D has 
been studied both along the continuous transition line, 
where critical exponents are computed, and in the coex- 
istence region, where FOPT line, latent heat curve and 



spinodal lines are reported. This latter observations con- 
firm the claim of Fernandez et al. [26j about the exis- 
tence of such transitions in quenched disordered short- 
range finite-dimensional systems. In the present model 
the FOPT can be seen by means of standard PT in the 
canonical ensemble, simply tuning an external pressure- 
like parameter (that is, a parameter adjustable in exper- 
iments on real samples [34]). Besides the peculiarity of 
FOPT, inverse freezing is also observed, for the first time, 
in a random short-range finite dimension system. Both 
features were absent in the same model on a hierarchical 



lattice 25 [ and in the ordered BC model no IT was ob- 
li. 



served 

Acknowledgments. Funding for this work has been 
provided by the INFM-CNR Seed Grant "Order in dis- 
order". We thank A. Nihat-Berker and H. Katzgraber 
for useful discussions. 



[1] J. Wilks and D.S. Betts, An introduction to liquid He- 
lium, Clarendon Press (Oxford, 1987). 

[2] C. Chevillard and M.A.V. Axelos, Colloid. Pol. Sci. 275, 
537 (1997). 

[3] S. Rastogi, G.W.H. Holme and A. Keller, Macro- 
molecules 32, 8897 (1999); A.L. Greer, Nature 404, 134 
(2000); N.J.L. van Ruth and S. Rastogi, Macromolecules 
37, 8191 (2004). 

[4] M. Plazanet et al. J. Chem. Phys. 125 154504 (2006). R. 
Angelini and G. Ruocco, Phil. Mag. B 87, 553 (2007). 
R. Angelini, G. Salvi and G. Ruocco, Phil. Mag. B 88, 
4109 (2008). R. Angelini, G. Ruocco, S. De Panfilis, Phys. 
Rev. E 78, 020502(R) (2008). E. Tombari et al, J. Chem. 
Phys. 123, 051104 (2005). 

[5] B. Donnio et al, Adv. Mater. 19, 3534 (2007). 

[6] N. Avraham et al, Nature 411 451 (2001). 

[7] G. Tammann, "Kristallisieren und Schmelzen", Metzger 
und Wittig, Leipzig (1903). 

[8] Generally speaking, at a fixed parameter externally tun- 
ing the interaction strength such as concentration, chem- 
ical potential or magnetic field. 

[9] P.E. Cladis, Phys. Rev. Lett. 35, 48 (1975): 39, 720 
(1977). H. Ozbek et al, Ph. Trans. 75, 301 (2002). A. 
Srivastava, D. Sa and S. Singh, Eur. Phys. J. E 22, 111 
(2007) and references therein. O. Portmann, A. Vater- 
laus, and D. Pescia, Nature 422, 701 (2003). 
[10] P.W. Arisz, H.J.J. Kauw and J.J. Boon, Carbohydr. Res. 
271, 1 (1995). 

[11] A. Haque and ER. Morris, Carbohyd. Pol. 22, 161 
(1993). 

[12] F. Sciortino, A. Giacometti and G. Pastore, Phys. Rev. 

Lett. 103 237801 (2009). 
[13] N. Schupper and N.M. Shnerb, Phys. Rev. Lett. 93 

(2004) 037202. 

[14] N. Schupper and N.M. Shnerb, Phys. Rev. E 72, 046107 

(2005) . 

[15] MR. Feeney., P.G. Debenedetti, and F.H. Stillinger, J. 

Chem. Phys. 119 4582 (2003). S. Prestipino, Phys. Rev. 

E 75, 011107 (2007). 
[16] A. Crisanti and L. Leuzzi, Phys. Rev. Lett. 95, 087201 



(2005). 

[17] M. Sellitto, Phys. Rev. B 73 180202(R) (2006). 

[18] S. K. Ghatak, D. Sherrington, J. Phys. C: Solid State 
Phys. 10, 3149 (1977). 

[19] A. Crisanti and L. Leuzzi, Phys. Rev. Lett. 89 (2002) 
237204; Phys. Rev. B 70 (2004) 014409. 

[20] The transition is first order in the thermodynamic sense 
and is not related to the so-called random first order tran- 
sition occurring in mean-field models for structral glasses. 

[21] L. Leuzzi, Phil. Mag. 87, 543-551 (2006). 

[22] M. Blume, V.J. Emery and R.B. Griffiths, Phys. Rev. A 
4, 1071 (1971). 

[23] H. W. Capel, Physica 32, 966 (1966); M. Blume, Phys. 

Rev. 141, 517 (1966). 
[24] D.M. Saul, M. Wort is and D. Stauffer, Phys. Rev. B 9, 

4964 (1974). A. Nihat Berker and M. Wortis, Phys. Rev. 

B 14, 4946 (1976). A. K. Jain and D. P. Landau, Phys. 

Rev. B 22, 445 (1980). 
[25] V.O. Ozcelik and A. N. Berker, Phys. Rev. E 78, 031104 

(2008). 

[26] L.A. Fernandez et al, Phys. Rev. Lett. 100, 057201 
(2008). 

[27] S. Caracciolo et al, Nucl. Phys. B403, 475 (1993). 
[28] H. G. Ballesteros et al, Phys. Rev. B 62 (2000) 14237. 
[29] M. Palassini, S. Caracciolo, Phys. Rev. Lett. 82, 5128 
(1999). 

[30] T.L. Hill, Thermodynamics of Small Systems, Dover 
(2002). 

[31] T. Jorg, Phys. Rev. B 73, 224431 (2006). H.G Katz- 
graber, M. Korner and A. P. Young, Phys. Rev. B 73, 
224432 (2006). M. Hasenbusch et al, Phys. Rev. B 76, 
094402 (2007); ibid. 76 184202 (2007). 

[32] E.K. Riedel and F.J. Wegner, Phys. Rev. Lett. 29, 349 
(1972). J. Zinn- Justin, Quantum Field Theory and Criti- 
cal Phenomena, Oxford University Press (Oxford, 1989). 

[33] Though Pl(p) is not exactly zero between the two peaks, 
cf. inset of Fig. [3] the determination of D C (T) with the 
present precision is robust against changes of po in a rea- 
sonably wide domain of p. 

[34] Although in numerical simulations changing the pressure 



5 



or, e.g., the bond dilution [2(|, or the relative probabili- 
ties of the random bond values [35| is technically equiv- 
alent, the latter two cannot be tuned in a single real 
thermodynamic experiment and require the preparation 



of several samples with different microscopic properties. 
[35] F.P. Toldin, A. Pelissetto and E. Vicari, J. Stat. Phys. 
135, 1039 (2009). 



